library(sasld)

# -------------------------------------------------------------------
# ATTENZIONE: la funzione sample e' stata modificata 
# a partire da R 3.6.0 per cui i "vecchi risultati"
# non posso essere riprodotti di default
# http://hostingwin.unitn.it/micciolo/PeM/01-rng.pdf
# -------------------------------------------------------------------
# se si vogliono riprodurre i risultati dei Laboratori
# con le versioni pił recenti di R va eseguita la seguente istruzione
# -------------------------------------------------------------------
RNGkind(sample.kind="Rounding")
# -------------------------------------------------------------------


# 7.1 ---------------------------------------------------------------
simDC.prop(0.54,4000,nrep=1)

simDC.prop(0.54,4000,nrep=10)

out <- simDC.prop(0.54,4000,nrep=100000)
summary(out)

sqrt(0.54*0.46/4000)

quantile(out,probs=c(0.025,0.975))

hist(out, xlab="proporzione campionaria",ylab="frequenza assoluta",main="Distribuzione campionaria della proporzione campionaria")
box()

plot(out)
# -------------------------------------------------------------------




# 7.2 ---------------------------------------------------------------
set.seed(123456)
simDC.mean(25,8.22,10,nrep=10)

set.seed(123456)
out <- simDC.mean(25,8.22,10,nrep=100000)
summary(out)

8.22/sqrt(10)
quantile(out,probs=c(0.025,0.975))

25 - quantile(out,probs=c(0.025,0.975))
1.96*8.22/sqrt(10)
qnorm(c(0.025,0.975),25,8.22/sqrt(10))




set.seed(123456)
out <- simDC.mean(25,8.22,100,nrep=100000)
summary(out)

2.59179/0.81961
sqrt(10)
quantile(out,probs=c(0.025,0.975))

25 - quantile(out,probs=c(0.025,0.975))
1.96*8.22/sqrt(100)
qnorm(c(0.025,0.975),25,8.22/sqrt(100))
# -------------------------------------------------------------------




# 7.3 ---------------------------------------------------------------
(40-20)^2/12
sqrt((40-20)^2/12)

set.seed(123456)
out <- simDC.mean(n=2,nrep=100000,family="unif",min=20,max=40)
summary(out)

s2 <- (40-20)^2/12; sqrt(s2/2)
quantile(out,probs=c(0.025,0.975))

qnorm(c(0.025,0.975),30,sqrt((400/12)/2))
qnorm(c(0.025,0.975),30,sqrt(s2/2)) # equivalente alla precedente

plot(out)
hist(out, xlab="media campionaria",ylab="frequenza assoluta",main="Distribuzione triangolare")
box()


set.seed(123456)
out <- simDC.mean(n=10,nrep=10000,family="unif",min=20,max=40)
plot(out)


#
# figura pag. 603
x <- seq(0,10,0.001)
plot(c(0,5),c(0,2), type="n", xlab="", ylab="", main="Distribuzione esponenziale")
r <- 1/2; y <- r * exp(-r*x); lines(x,y,lty=3,lwd=2)
r <- 1;   y <- r * exp(-r*x); lines(x,y,lty=2,lwd=2)
r <- 2;   y <- r * exp(-r*x); lines(x,y,lty=1,lwd=2)
legend(1,2,c("rate = 2","rate = 1","rate = 0.5"),lty=c(1,2,3))


set.seed(123456)
out <- simDC.mean(n=10,nrep=10000,family="exp",rate=2)
plot(out)

set.seed(123456)
out <- simDC.mean(n=1000,nrep=10000,family="exp",rate=2)
plot(out)


# dati "reali" (distribuzione empirica)
pop <- c(1,2,2,2,3,3,3,3,3,4,4,4,
         4,4,4,4,5,5,5,5,5,5,5,5,
         5,6,6,6,6,6,6,6,6,6,6,6)
mean(pop)
s2 <- var(pop)*35/36
s <- sqrt(s2)

set.seed(123456)
out <- simDC.mean(n=30,nrep=10000,family="real",data=pop)
summary(out)
plot(out)

set.seed(123456)
out <- simDC.mean(n=300,nrep=10000,family="real",data=pop)
summary(out)
plot(out)